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Abstract 

We introduce a one-parameter family, < H < 1, of pair potential functions with a single relative 
energy minimum that stabilize a range of vacancy-riddled crystals as ground states. The "quintic 
potential" is a short-ranged, nonnegative pair potential with a single local minimum of height H 
at unit distance and vanishes cubically at a distance of y/3. We have developed this potential 
to produce ground states with the symmetry of the triangular lattice while favoring the presence 
of vacancies. After an exhaustive search using various optimization and simulation methods, we 
believe that we have determined the ground states for all pressures, densities, and < H < 1. For 
specific areas below 3^/3/2, the ground states of the "quintic potential" include high-density and 
low-density triangular lattices, kagome and honeycomb crystals, and stripes. We find that these 
ground states are mechanically stable but are difficult to self-assemble in computer simulations 
without defects. For specific areas above 3\/3/2, these systems have a ground-state phase diagram 
that corresponds to hard disks with radius For the special case of H = 0, a broad range 
of ground states is available. Analysis of this case suggests that among many ground states, a 
high-density triangular lattice, low-density triangular lattice, and striped phases have the highest 
entropy for certain densities. The simplicity of this potential makes it an attractive candidate for 
experimental realization with application to the development of novel colloidal crystals or photonic 
materials. 

PACS numbers: 
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I. INTRODUCTION 



The field of self-assembly provides numerous examples of using atoms, molecules, colloids, 
and polymers as building blocks in the development of self-organizing functional materials. 
For example, researchers have used nanoparticles to create stacked rings and spirals us- 
ing DNA linkers, 1 to construct photonic crystals using binary systems of colloids,- and to 
prototype a self-organzing colloidal battery.— 

Recently, there have been significant efforts to design pair potentials that yield self- 
assembly of targeted ground-state (potential energy minimizing) structures using "inverse 
methods. With an inverse approach, one chooses a targeted structure or property and 
uses optimization methods to develop a robust pair potential function. We envision that 
colloids will be the experimental testbed for the optimized interactions because colloids 
offer a significant range of repulsive and attractive interactions that can be tailored in the 
laboratory- For these models, the many-body interactions of a system are reduced to pair 
interactions. For a pair potential v(r), where r is the distance between two particles, the 
potential energy per particle u of a many-body configuration r N of N particles reduces to 



Inverse approaches have been applied to produce low-coordinated ground states such as 
honeycomb, 7 square, 8 simple cubic,' diamond and wurtzite crystals^ as well as many-body 
configurations on a sphere. 5 In addition to targeting material structures, the inverse approach 
has been applied to material properties including negative Poisson's ratio,— negative thermal 
expansion,— and scattering characteristics of many-particle systems.— The ability to control 
vacancy concentrations and arrangements in ground-state structures is one property yet 
to be achieved via inverse methods. We expect the design of vacancy-riddled crystals to 
have important technological applications and may provide fundamental insight into certain 
physical phenomena. The scattering of light, 1 - ionic conductivity, 15 transport processes in 
heterogeneous materials,- and possibly supersolid behavior in quantum systems^ 7 - 1 ^ are 
directly related to vacancies in a material structure. 

In this paper, we introduce the "quintic potential" as a pair interaction function that 
yields vacancy-riddled lattices and low-coordinated crystals as ground states at high density 
and striped patterns as ground states at low-density. Although not a result of a true inverse 
method, this family of potential functions arose through a selection process so that vacancy- 
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FIG. 1: (Color online) Quintic potential, Eq. [6] for various H values. The scales of energy and 
length are dimensionless. 

riddled and low-coordinated lattices are energetically degenerate or possibly favorable to the 
triangular lattice for certain densities. Each quintic potential function curve, illustrated in 
Fig. [TJ is steeply repulsive for small r, has a local minimum at r = 1 so that v(l) = H, is 
nonnegative, and is smoothly truncated to be zero at r = \/3 and beyond. The algebraic 
form of the potential is given in Sec. [TTJ The understanding of the phase behavior associated 
with the quintic potential is the first step toward developing a more comprehensive inverse 
approach to controlling vacancies in ground-state structures. 

We find that ground states can vary from a high-density triangular lattice, the kagome 
and honeycomb crystals, a striped phase, a low-density triangular lattice, and a hard-disk- 
like fluid as the pressure decreases. The positivity of the local relative minimum is a key 
feature that allows vacancies to be stabilized in the ground state. Consider a system of 
particles arranged in a triangular lattice at a specific area (area per particle) a = V3/2 so 
that each of the particle's first neighbors lie at r = 1 and second neighbors at r = \/3. 
Because the local minimum in v (r) is at positive energy, the removal of a particle will lower 
the total potential energy of the system. However, to maintain the constraint on a, the 
lattice spacing must decrease, which increases the total potential energy. This family of 
potentials is designed so that the creation of a vacancy and the subsequent reduction in 
lattice spacing yields a net decrease in the potential energy per particle. 
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However, the general control of vacancies is a challenging task due to long-ranged vacancy- 
vacancy interactions. For the Lennard- Jones and many short-ranged repulsive potentials, 
vacancies arise in equilibrium solids at positive temperature due to their entropic contribu- 
tion to the free energy— At T = 0, they can be "frozen in" during cooling, though they 
are not present in the equilibrium ground states of these systems. For typical potentials 
whose ground states are ordered, some pair potentials, such as the Lennard- Jones interac- 
tion, vacancies cost energy due to missing "bonds" of negative energy, while with repulsive 
potentials and the electron crystal, vacancies cost energy due to strain energy in the crystal. 

If a vacancy is introduced into an otherwise perfect crystal at positive pressure, strains 
will arise so that the system can reduce its potential energy and achieve mechanical stabil- 
ity. In typical materials such as Lennard- Jones systems, these relaxations cannot reduce the 
energy below that of the ground state. When more than one vacancy is present, the strain 
fields interact and subsequently mediate long-ranged vacancy-vacancy interactions. These 
interactions are highly nontrivial even in the case of two vacancies in a crystal. In two 
dimensions, several studies have examined these effective interaction energies between two 
vacancies as a function of separation distance for a number of potentials including the elec- 
tron crystal (Coulombic potential with positive background charge),— - — Lennard- Jones,— 
Gaussian core, 24 and a repulsive r -3 potential.— 

In the two-dimensional electron crystal, vacancy-vacancy interactions are attractive and 
fall off as r~ 3 ,— where r v is the separation distance between two vacancies, which agrees 
with continuum elasticity theory— For the Lennard- Jones potential, vacancy-vacancy inter- 
actions are also attractive at least for short distances.— These effective vacancy interaction 
potentials are not necessarily monotonic as a function of lattice spacing.—^ The displace- 
ment fields can be highly anisotropic near a vacancy.—^ In the vicinity of the vacancy, the 
atomistic details account for the discrepancy between numerical results for the displacement 
fields and those predicted by continuum elasticity theory. Beyond fifteen lattice spacings, 
this continuum theory accurately reproduced the results of numerical studies under some 
boundary conditions.— 

The arrangement of low concentrations of vacancies may not simply be controllable be- 
cause of the nontrivial coupling between the pair potential function, the local deformations 
that it produces when a vacancy is introduced into a perfect crystal, and the effective 
vacancy-vacancy interactions. However, a number of pair potential functions have been 
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found to give rise to ground states with high concentrations of vacancies. For example, the 
"honeycomb potential" was developed via inverse methods to yield the honeycomb crys- 
tal as a robust ground-state structure.-^ At positive temperatures, the honeycomb struc- 
ture remained stable for a large temperature and pressure region of the phase diagram.— 
This double-well potential was generalized as the Lennard- Jones- Gaussian potential, which 
showed a number of complex crystal structures including pentagonal, hexagonal, nonago- 
nal, and decagonal unit cells as the depths and locations of the energy wells were varied.—^ 
Colloidal particles with three "patches," or attractive interaction sites on the surface of a par- 
ticle, have recently been shown to yield a honeycomb structure for some pressures at T = 
in addition to other low-coordinated crystal structures.— In addition, the square-shoulder 
potential, a hard-core potential with a soft corona of variable length, shares some of the same 
ground-state features as the quintic potential. In particular, low-coordinated structures such 
as a honeycomb-like trivalent configuration and lanes are ground-state features 3 - 1 " 1 ^ that are 
shared with the quintic potential. Such a potential can also form lamellar and micellar 
phases for relatively long-ranged coronas as was shown theoretically.— 

Lastly, we note that several potentials qualitatively similar to the quintic potential have 
been studied in other contexts. For example, a potential that gives rise to quasiperiodic 
order and the square lattice in two dimensions has one local minimum and local maximum 
and finite cutoff, 34 as does the Dzugutov potential 3 ^ which gives rise to a number of com- 
plex local clusters in three dimensions.— However, both of these potentials have a local 
attractive minimum (i.e. negative potential energy). A piecewise linear potential with a 
hard core, single local minimum, and single local maximum showed the formation of chains 
and labyrinths at positive temperature, though the ground state is not fully characterized.— 
However, in contrast to the above potentials, the quintic potential is short-ranged, positive, 
isotropic, continuous and differentiable, which may be beneficial to experimentalists hoping 
to realize this potential in a laboratory setting. We find that with this potential, a number 
of low-coordinated structures arise as ground states including the honeycomb and kagome 
crystals as well as a "striped" phase. This represents a significant achievement and a first 
step toward controlling vacancies in ground-state structures. 

The remainder of this paper is as follows. The quintic potential is defined and detailed in 
Sec. [TT]while our methodology is detailed in Sec. IIHI In Sec. HV] we define the phase diagram 
in the specific area- if plane and in the pressure- if plane. We discuss the characteristics of 
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the phases, the results of simulation, and their mechanical stability. In Sec. |Vj we focus on 
the special case where H = 0. This case yields anomalous behavior because the potential 
energy vanishes inside and outside the energy barrier. Lastly, we discuss the implications of 
this potential and identify extensions to this work in Sec. I VI I 



II. QUINTIC POTENTIAL 



The generalized quintic potential consists of a fifth-order polynomial that is truncated 
beyond \/3 to be zero. It has the form 

5 / _\ 4 / x 3 



f(r) 



L{m) (r — Vtj + K(m) (r-Vz) - (r - V?> 



< V3 



(1) 



and zero otherwise, where m is an unsealed height of the local minimum. The coefficients 
L(m) and K(m) are chosen so that r = 1 is a local minimum and f(l) = m, and, as a 
function of m, are given respectively as 

-5 / _ \ -2 

(2) 
(3) 



L(m) = Am (Vz - ij ° - (V3 - 1 
K{m) = 5m(V3-l) ^ - 2 (v^ - l) . 



The condition that /"(l) > ensures that the stationary point is a relative minimum, and 
therefore the parameter m must be constrainted to 

(4) 
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The position of the energy barrier (the relative maximum) occurs at 
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The generalized pair potential is rescaled so that u(r&) = 1 to set a uniform energy scale, 

v(r) = f(r)/f(r b ). (6) 

Upon rescaling, the height of the relative minimum H is defined so that v(l) = H and 
varies from zero to unity. When H — 1, the potential is monotonically decreasing and is flat 
at r = 1. In this case, there is no local minimum or maximum. These soft potentials are 
repulsive for small r, and the first and second derivatives vanish at r = a/3. We construct 
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FIG. 2: The relationship between m and H for the quintic potential. 

the potential such that the first and second derivatives vanish at r = \/3 so that the second- 
order expansions of the potential energy required for stability calculations {e.g. phonons) are 
well-defined. Examples from the family of quintic potentials are shown in Fig. [TJ Because 
developing a simple expression relating m and H is difficult, a table relating m and H was 
used. The relation between m and H is shown in Fig. [2J 

III. METHODS 

At T = 0, the free energy per particle g is identical to the enthalpy per particle and 
is related to the average potential energy per particle u via g = u + pa. To map the 
ground-state phase diagram as a function of specific area a, pressure p, and H, we utilized 
the double-tangent construction method to find the lowest free-energy structure among 
candidate ground-state configurations. The double-tangent construction is demonstrated in 
Fig. El 

We considered several crystal structures as candidate ground states including the trian- 
gular lattice (TRI), square lattice (SQ), honeycomb crystal (HC), kagome crystal (KAG), 
and its close relative, the "anti-kagome" crystal (AKG). Whereas the kagome crystal has 
vacancies located on a triangular lattice separated by two nearest neighbor spacings, the 
anti-kagome crystal has vacancies located on a rectangular lattice. For all densities, the 
anti-kagome has a potential energy greater than or equal to the kagome lattice due to the 
differing local coordination numbers. We considered these crystal structures because they 
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are relative simple crystal structures with varying degrees of local coordination. The con- 
struction of the potential, with a well at r = 1 and vanishing energy at r = \f?> would 
intuitively favor such geometries. We use other techniques to systematically explore other 
crystals with a small number of particles in the unit cell. With these crystal structures, we 
performed lattice sums over the relevant range of specific area. 

We also performed an exhaustive search for energy-minimizing n-particle crystals, peri- 
odic configurations containing n particles per unit cell, while allowing for shape deformation 
of the unit cell. For example, consider a two-particle crystal with lattice vectors [d 1: 0] and 
[<i 2 , d 3 ] and a basis where particle is at (0, 0) and particle 1 is located at (xi,yi). For this 
system, the potential energy per particle is given as 

u = v(r 01 ) + K r oo') + v(r iv ) + u(r„i') + v(t 10 >)] , (7) 

i 

where the summation is over all lattice sites except the origin site and 0' represents particle 
in a cell other than the origin cell (and likewise for particle 1). Minimizing u for a fixed 
area per particle a and eliminating d 3 = na/di casts this as an unconstrained minimization 
that is function of d±, d 2 , X\ and y\. The function is minimized using the conjugate gradient 
algorithm. It is straightforward to generalize this to a larger number of particles in the unit 
cell. 

To ensure that we obtain globally optimal solutions, we use a large number of initial 
conditions slowly and systematically varying d 1 , d 2 , X\ and y\ from zero to ny/3. After each 
minimization, we ensure that the unit cell is not significantly sheared. A highly sheared 
unit cell requires summation over a large number of unit cells to ensure that all interactions 
are accounted for and therefore we discard those solutions with angles less than 10°. This 
procedure is repeated for nearly two thousand values of a on the range 0.7 < a < 3^3/2 so 
that a smooth u versus a curve is generated. 

These minimizations are a function of In variables, and therefore many minimizations 
at a fixed area per particle can be performed easily. However, as the number of particles 
per unit cell grows, the number of initial conditions required to ensure global optimality 
grows as 2 n . We performed these minimizations systematically for up to four particles per 
unit cell, beyond which it becomes too intense to systematically explore thousands of initial 
conditions for thousands of specific areas. 

We also used slow cooling via molecular dynamics (MD) to obtain candidate ground-state 
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structures for larger and possibly disordered systems. With a system of 800 particles in a pe- 
riodic box, the system was simulated using the Verlet algorithm and Andersen thermostat.— 
The system was initialized as a high-temperature liquid and slowly cooled according to a lin- 
ear temperature schedule from T = 0.4 to 0.025 over the course of 1.6 xlO 7 time steps. After 
the simulation was terminated, the system was quenched to a potential energy minimum 
using the conjugate gradient algorithm. 

Lastly, we performed lattice Monte Carlo (MC) optimizations using several hundred 
lattice sites. In these optimizations, the specific area was fixed and the lattice sites were 
either occupied or unoccupied by particles. The optimization variables were the occupation 
state of the lattice sites, the angle between the lattice vectors and the length of one lattice 
vector. The length of the other lattice vector was constrained by the area, number of 
particles, angle between the lattice vectors, and the length of the first lattice vector. Possible 
Monte Carlo moves included inserting a particle to a vacant site, removing a particle from 
an occupied site, swapping a particle from an occupied site to an unoccupied site, perturbing 
the angle of the lattice, and perturbing the length of one of the lattice vectors. The system 
was then optimized via simulated annealing. It was given an initial "temperature," or energy 
scale, and moves were accepted or reject based upon the Metropolis algorithm. Ten thousand 
MC trials were attempted in each cycle and several hundred cycles were performed for each 
optimization. 

It is important to note that although many interesting structural characteristics arose 
from the results of molecular dynamics simulation and Monte Carlo optimization, no fi- 
nal configurations were identified as ground-state structures. The lattice sums and crystal 
optimization methods always provided the lowest free-energy structures. 

The double-tangent construction requires care and precision since other phases can come 
very close to the coexistence free energy. The p-g and u-a curves are discretized over several 
thousand data points and linearized between between data points. The identification of 
coexisting phases can be done by using the lower, concave envelope of the p-g curve. The 
coexistence points and ranges of stability were then calculated using a tabulation of u, a, p 
and g and linear interpolation. 
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FIG. 3: (Color online) Potential energy per particle u as a function of specific area a for selected 
crystals and the optimal 2- and 3-particle crystals for H = 0.5. The dense triangular, honeycomb, 
striped, and open triangular crystals are stable ground states for H = 0.5. The square and anti- 
kagome structures are omitted for clarity. Dashed lines represent the double tangent construction. 

IV. RESULTS 



A. Phase Diagram 

The double-tangent construction revealed two different types of phase behavior depending 
on the value of H. In no case did the configurations resulting from MD or MC yield a ground 
state. For < H < 0.762902, there exist five distinct structures. Figure |3] illustrates the u-a 
curves for H = 0.5. In the figure, the optimal n-particle crystals are not visible when they 
are identical to the triangular, honeycomb, or kagome crystals. As shown in in Fig. [31 at 
highest pressure, a dense triangular lattice (TRI), where neighboring particles lie inside the 
potential energy barrier of other neighboring particles, is the ground state. At a reduction 
of pressure to p — 1.84432, the dense triangular lattice is in coexistence with the honeycomb 
crystal (HC), which is illustrated with the dashed tangent lines. In Fig. El the curve for the 
kagome lattice appears possibly to touch the coexistence line due to the finite thickness of 
the line. However, plotting p against g reveals an intersection between the curves for the 
TRI and HC structures. The kagome lattice has a higher free energy than both structures 
in this density range. 
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Further reduction of the pressure to 0.75263 yields a "striped" (ST) or lane phase coex- 
isting with the honeycomb phase. This phase is an affinely-stretched triangular lattice. The 
first few coordination shells contain two, four, and two particles. The curves representing 
the lowest-energy crystals with a two- and three-particle basis either the triangular lattice, 
honeycomb crystal, or kagome crystals generally for a < 1.3 with a corresponding basis. 
For larger a, these curves appear to come close to the coexistence line between the HC and 
ST phases or the ST and open TRI phases. The p-g curves reveal that the free energy is 
always above that of the coexisting phases. These near-ground states are generally crystals 
that nearly mimic the coexistence. For example, the optimal three-particle crystal between 
the HC and ST phases is an alternating sequence of lines with honeycomb-like coordination 
and striped coordination. These are suboptimal structures because the specific neighbor 
lengths do not match exactly. At sufficiently low pressure, for H = 0.50 and p = 0.58275, 
the ST phase coexists with an "open" triangular lattice where the neighbor particles lie on 
the outside of the potential energy barrier. 

Lastly, at vanishing pressure, any configuration in which particles are separated by at 
least a/3 is a ground state since the free-energy vanishes. This occurs for a > 3^3/2. In 
this region, the ground states behave like hard disks (HD) of radius a/3. There would be 
a hard-disk crystal and liquid regime with specific area scaled to the hard-disk equation of 
state, which has been well characterized.^, 

For 0.762902 > H > 1, the kagome crystal emerges as a ground state structure for a nar- 
row range of stability. This is an especially important result. In previous work, researchers 
used an inverse methodology to design pair interactions for targeted ground states.- Al- 
though they were successful in engineering potential for the honeycomb and square crystals, 
they were unable to do so for the kagome. The area and pressure ranges of stability increase 
with H. The emergence of the kagome crystal as a ground state is shown in the p-a curve 
for H = 1.00 in Fig. HJ As the pressure is reduced, the sequence of stable phases is the dense 
TRI, KAG, HC, ST, open TRI, and HD phases. 

The full phase diagrams in the p-H plane and the a-h plane are shown in Figs. and [61 
respectively. Figure [5] shows that the coexistence lines are nearly linear for small H . The 
pressure range of stability for the open TRI, ST, and KAG phases widens as H increases. 
When H is sufficiently large to include the KAG phase, the pressure range of stability for 
the HC phase narrows. As H increases, the area range of stability increases for most phases, 
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FIG. 4: (Color online) Potential energy per particle u as a function of specific area a for selected 
crystals and the optimal 2- and 3-particle crystals for H = 1.00. The dense triangular, kagome, 
honeycomb, striped, and open triangular crystals are stable ground states for H = 0.5 while for 
H = 1.00, the kagome also emerges as a stable ground state. The square and anti-kagome structures 
are omitted for clarity. Dashed lines represent the double tangent construction. 

excluding the dense TRI phase, which is evident in Fig. [6j 

Near H = 1.00, the slopes of the curves change dramatically. This is due to the softening 
of the pair potential function as H approaches unity and the fact that the relative minimum 
and maximum come together much more rapidly as H approaches unity. Fig. [1] shows that 
the H = 1.00 potential is significantly softer near r = 1 than other potentials. 

The softening of the potential and the rise of the relative minimum are important features 
for the inclusion of the KAG phase. In comparing Fig. [3] to Fig. HI one can see that the curve 
for the kagome crystal initially begins above the double tangent connecting the TRI and HC 
phases. As the potential softens and the relative minimum increases, the first minimum in 
each curve increases proportionally to H. However, the softening of the potential bends the 
left-side of the curves in such a way that the KAG curve lies at lower free energy than that 
TRI-HC coexistence line. The combination of the height of the relative minimum and the 
softening are directly related to the potential energy and pressure of the system, the two 
components of the T = free energy. 

The KAG phase emerges at a triple point at H = 0.762909 and p = 2.938, where the 
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FIG. 5: (Color online) Phase diagram in the p-H plane as calculated by double tangent method. 
The hard-disk-like state (HD) occurs for p = and H > and the filled triangle on the lower curve 
represents the triple point. 
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FIG. 6: Phase diagram in the a-H plane as calculated by the double tangent method. 
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FIG. 7: (Color online) Comparison of the location and coordination number of the triangular 
lattice, honeycomb crystal and kagome crystal for a = 0.82443, 1.18745, and 1.06445, respectively, 
for H = 0.762902. At this H and p = 2.938, the TRI, HC, an KAG phases form a triple point and 
are in coexistence as ground-state structures. The delicate balance between potential energy and 
forces allow for coexistence. 

TRI, HON, and KAG phases form a triple point. We have examined the subtle interplay 
between the shape of the pair potential function and its derivative at this coexistence point 
in Fig. [71 The figure shows v (r) and —v'(r)/2 for H = 0.762902. The plot also shows the 
location of the nearest and next-nearest neighbors for the crystals. The HC phase has the 
closest neighbors while the TRI phase has the smallest forces. 

The coexistence of these structures requires that a complex system of equations be sat- 
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isfied. For the coexistence pressure p* and free energy g*, the system of equations becomes 



P 
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g* 

p* 



v(d t ) + v{d t V3) - ^v'(d t ) - ^p-v'(d t V3) 
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(8) 

(9) 
(10) 

(11) 
(12) 

(13) 



where di is the nearest-neighbor distance for each structure and is necessarily less than unity 
for the quintic potential. Evidently, the quintic potential for H = 0.762902 satisfies these 
equations. The complexity of the coexistence equations is clear in that each structure weighs 
certain parts differently. The design of other potentials that stabilize these structures may 
utilize this system of equations within an appropriate optimization framework. 



B. Stability 

The positivity of the squared frequency of all phonons ensures mechanical stability of a 
crystal. We have examined the phonon spectra for the KAG, HC, and ST phases. These 
phases are confirmed to be mechanically stable within the relevant density and H ranges. 
Figure |8] illustrates the phonon spectra for certain points in the reduced first Brillouin zone 
for the KAG, HC, and ST structures for H = 0.875. THE HC and KAG structures have one 
particularly soft acoustic branch, labeled in the figures. Using the ratio ^max/^min a ^ M 
point a simple metric for the extent of mechanical stability, this ratio stays relatively fixed 
for the KAG lattice as H increases. This marks an achievement since previous attempts 
to stabilize the kagome crystal had been unsuccessful.- For the HC structure, this ratio is 
higher for H = 1.00 than for that which is plotted in Fig. |8l indicating that the crystal is 
more stable as H increases. In general, the quintic potential and the honeycomb potential 8 
yield similar ratios. The striped phase whose lattice vectors for a = 1.51 and H = 0.875 are 
[1.65029, 0] and [0.27688, 0.91490] is also mechanically stable. The phonon spectra from the 
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(a) (b) (c) 

FIG. 8: (Color online) Phonon spectra for (a) the kagome crystal at a = 1.035 and (b) the 
honeycomb crystal at a = 1.2, and (c) the striped phase at a = 1.51 for wave vectors in the 
reduced first Brillouin zone for H = 0.875. The first Brillouin zone for the striped phase is shown 
in inset in (c). These crystals are mechanically stable at these densities with the quintic potential. 
The branches containing soft acoustic modes are denoted with a *. 



other "quadrants" are identical to the quadrant displayed in Fig. 8(c) due to the symmetry 
of the Brillouin zone. 



C. Low-Energy Structures 

Although those structure obtained by MD and MC methods were suboptimal, Figures [3] 
and |H show that the free-energy differences between these structures and the ground states 
are small. Several structural motifs emerge for this potential and are shown in Fig. for 
H = 0.5. These metastable states may have technological value if the freezing kinetics are 
sufficiently slow. Preliminary simulations suggest the freezing behavior varies with density. 

For example, with a < 1.2, systems typically freeze into a triangular lattice with vacancies 



randomly distributed as in Fig. 9(a) Although this is not a ground-state structure, it yields 
a vacancy-riddled lattice. The vacancies appear to have no particular order. The freezing 
transition, the temperature at which the system changes from a high-density liquid to an 
ordered phase, appears to be first order. The drop in potential energy as the temperature 
is reduced is sharp in this density range. 

Systems at densities where coexistence between the HC and ST phases are the ground 
states tend to exhibit a freezing from the liquid phase to a rigid, structured phase. However, 
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the structured phase, which we believe to be metastable, is characterized by rings and strings 



as in Fig. 9(b) The six-particle ring is a characteristic of the HC structure and the strings 
are characteristics of the ST phase. These metastable states are disordered due to the fluid 
nature of the strings. The drop in potential with temperature is not as sharp for these 
densities compared to those at higher densities. 

Lastly, at lower densities, as with those associated with the ST and open TRI phases, the 



metastable, low-energy states adopt labyrinthine characteristics, Fig. 9(c), and eventually 



colloidal polymers and monomers, 9(d) In this density range, the drop in potential energy 
associated with freezing is weak. The phenomena and characteristics of low-energy states 
are common to all H. However, as H increases so does the propensity to form labyrinthine 
characteristics. This is due to the lower energy barrier and the decrease in r&, the location 
of the energy barrier. A qualitatively similar, but piecewise linear, potential also yields 
a labyrinth phase at positive temperature.— Although the positive temperature behavior 
has not been fully explored in this paper, we anticipate that the equation of state for this 
family of potentials will have interesting behavior. Manipulation of the kinetics to achieve 
such unusual structures represents one path to achieving kinetically stable materials with 
controlled vacancy concentrations. 



V. SPECIAL CASES: H = 



The H = pair potential has interactions that vanish at a pair distance of unity and 
beyond y/3. At positive pressure, the ground state is the dense triangular lattice. However, 
at zero pressure and a > v3/2, ground-state configurations will have a vanishing poten- 
tial energy and pressure (enthalpy vanishes). Thus, the type of available ground states is 
dependent on the area. 

A number of interesting structures arise as ground states. For a > v3/2, dilutions of the 
triangular lattice with unit neighbor spacing are ground states, including the honeycomb 
and kagome crystals and other lattice gases. For a = vTl/2, the striped phase, a Bravais 
lattice with lattice vectors of [1,0] and [1/2, \/Tl/2], becomes available as a ground state. 
Each particle has two neighbors at unit separation and four neighbors at separation of \/3. 
For a > ^/Tl/2, dilutions of this lattice are also ground states. In addition, any small 
expansion in the direction normal to the stripes, will allow each stripe to gain some fluidity. 
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(c) a = 1.7 (d) a = 2.3 



FIG. 9: (Color online) Metastable configurations obtained via molecular dynamics for H = 0.5. 
In an attempt to cool to a liquid to the ground state, the system has difficulty crystallizing 
and/or phase separating into the appropriate ground-state structures. As density is decreased, 
the metastable states go from ordered structures with vacancies to a labyrinthine network to a 
monomer ic liquid. 

Complex liquids of strings have found as ground states using molecular dynamics for at least 
a > 2.15. However, the geometric problem is highly nontrivial since for several runs with 
a > 2.15, ground states could not always be obtained. For a = 3v3/2, an open triangular 
lattice with neighbor spacing of V3 is a ground state structure. 

The question remains as to which type of configurations are entropically (thermodynami- 
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cally) favorable, or rather which type of system makes up the largest fraction of configuration 
space. We first make the distinction between discrete (e.g. lattice gas) and continuous en- 
tropy (e.g. hard-disk crystal). In the classical case, continuous entropy is uncountable while 
discrete entropy is not. Making the analogy to hard disks and hard spheres, we believe that 
the highest entropy phases for a specified area are those that have the availability for contin- 
uous entropy. Configurations with the highest-dimensional configuration space dominate the 
T = entropy. We estimate that, for a certain a, configurations with the fewest constrained 
degrees of freedom will have the highest dimensional configuration space, and therefore the 
highest entropy. For a < y/3/2, all degrees of freedom per particle are constrained. For 
a > 3v3/2, no degrees of freedom per particle are constrained. 

In order to gain continuous entropy, the system must have the ability to expand the 
a/3 "bonds," which can be either first-neighbor or second-neighbor connections between 
particles. Any slight expansion of the dense triangular lattice is not a ground state because 
each neighbor must be constrained to unit separation. However, an expansion of the open 
triangular lattice is still a ground state because the v3 bonds are not constrained from 
above. The problem can be cast as a tiling problem of four triangular tiles. The possible 
triangular tiles are those with side lengths of unity or v3. These are depicted in Table |T] 
along with their shorthand notations. 

Because there are no constraints beyond distances of a/3, the edges with lengths of a/3 
can be considered "elastic." For simplicity, we first consider the tiling problem where these 
edge lengths are fixed. The rules of the tiling problem are as follows: 

• The plane must be tiled with no gaps. 

• For pairs of tiles that share an edge, the edges must be of the same length so that the 
vertices also match. 

• The (l,l,r) and (l,r,r) tiles cannot share a a/3 edge since this violates a second neighbor 
constraint (i.e. two vertices are separated by a length between unity and a/3). 

The internal angle of the (l,r,r) tiles are such that they cannot integrate well with the 
other tiles. Therefore, in any tiling that includes (l,r,r) tiles, these tiles must "phase sepa- 
rate" from the others. We consider two types of tilings - those with (l,r,r) tiles and those 
without. First, we consider those without (l,r,r) tiles. An example of such a tiling is shown 
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TABLE I: Available triangular tiles for the H = tiling problem. 



Name 


Shape 


Angles (°) 


Area 




/\ 


fin fin fin 


\/3/4 


(r,r,r) 


/\ 


6n, 60, 60 


3^/4 


(l,l,r) 




30, 30, 120 


V3/4 


(l,r,r) 




33.56, 73.22, 73.22 


vTI/4 



in Fig. 10(a) where each vertex is decorated with a particle. These tilings are simply di- 
lutions of the dense triangular lattice, or coexistence between a dense and open triangular 
lattice. The (1,1, r) tiles have two roles. They can act as intermediaries between domains 
of the dense triangular lattice and the open triangular lattice, and they can dimerize to 
form a rhombus making a domain of the dense triangular lattice. This coexistence between 
the dense and open triangular lattices is represented as a double-tangent line connecting 
the fully constrained dense triangular lattice to the unconstrained open triangular lattice as 
shown in Fig. [TTJ 

Next we consider those tilings that include (l,r,r) tiles. Two periodic tilings with a 



specific area of a = yll/2 « 1.658312 exist and are shown in Figs. 10(b) and 10(c), We 
deem these the stripe and "zig-zag" tilings respectively. In these phases, any small expansion 
perpendicular to the stripe or zig-zag will allow the stripes to have fluidity. On average, 
these configurations will constrain one degree of freedom per particle, as plotted in Fig. [TTJ 
For any tiling consisting of all tiles, the (l,r,r) tiles must segregate so that the the other 
tiles can fully tile the plane. Dimerizing the (l,r,r) tiles along the \/3 edges creates a wide 
stripe. Building off the wide stripe requires (1,1,1) tiles or (l,l,r) tiles and forms a local 



coexistence with the dense triangular lattice as shown in Fig. 10(d) Dimerizing the (l,r,r) 
tiles along the unit length edge creates a narrow stripe. The exposed edges have length \/3, 
requiring (r,r,r) tiles to be directly adjacent to the stripe. This forms a coexistence between 
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(a) (b) (c) 




(d) (e) 

FIG. 10: (Color online) (a) Tiling excluding (l,r,r) tiles, (b) striped phase using only (l,r,r) tiles, 
(c) zig-zag phase using only (l,r,r) tiles (d) coexistence between dense triangular lattice and the 
striped phase, and (e) coexistence between the open triangular lattice and the striped phase. 

the the striped phase and the open triangular lattice. The double tangent construction of 
these coexistences is displayed as the solid red line in Fig. CLU These coexistences between the 
striped phase and the dense triangular lattice and the striped phase and the open triangular 
lattice maximize the number of unconstrained degrees of freedom and likely maximize the 
entropy of the ground state. In addition, these system can take on discrete entropy by way 
of stacking variants of striped phase and the lattice phase. Using "S" and "T" to denote a 
stripe and a triangular lattice layer, a STTTTSSS system is degenerate with a STSTSTST 
system. 

The zig-zag phase can only form local coexistence with the open triangular lattice. To do 
so, the (l,r,r) tiles must form a trimer which exposes the a/3 edges on either side of the trimer. 
Since they are in trimers, their stacking entropy would be less that the stacking entropy 
available to the coexistence between the stripe and open triangular lattice. Therefore, we 
believe the stripe phase has a higher entropy than the zig-zag phase. 
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a = A/N 

FIG. 11: (Color online) Unconstrained degrees of freedom per particle for candidate ground-state 
structures for H = potential. Configurations with the maximum number of unconstrained degrees 
of freedom are presumed to be the ground state (red circles). 

VI. DISCUSSION AND CONCLUDING REMARKS 

In this paper, we have developed the ground-state phase diagram of a new potential 
that gives rise to a number of novel phases that include low-coordinated crystal structures. 
Given the unusual nature of the ground state, we expect that the equation of state and full 
phase diagram for these systems to exhibit other unusual behaviors at positive temperature. 
For example, in the global phase diagram of the the Lennard- Jones-Gauss, or honeycomb 
potential, there was no gas-liquid coexistence in the low-temperature, low-density part of 
the phase diagram, nor was there a liquid-liquid phase coexistence^. We expect most of 
the solid-solid transitions that we find at zero temperature will remain at small nonzero 
temperature. Our preliminary calculations suggest that strings, or polymers, may arise in 
equilibrium at low densities. 

In addition, we have interest in the mobility of vacancies, particularly in the cases where 
vacancy concentrations are dilute, (e.g. H = and a m a/3/2), due to the possible relation 
to supersolid behavior.—^ For a just above a/3/2 and H = 0, we expect there to be a small 
number of vacancies in the system at low-temperature. We have estimated the potential 
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energy required for a particle to "hop" from a lattice site to a vacant site. By initializing a 
system with one vacancy with a "hopping" particle along the transition to the vacant site, 
in an otherwise undistorted lattice, we minimized F(r N ) = |Vu(r Ar )| 2 , where represents 
particle coordinates, using conjugate gradient minimization. The resulting configuration 
represents a saddle point in the energy landscape. The difference in potential energy of the 
saddle point and the ground state is the energy barrier required for the particle and vacancy 
to swap positions. This barrier sets an activation energy for classical thermal motion. For 
the quantum case, it would also enter in the tunneling rate. For H = potential, there is a 
single saddle point in the vicinity of the numerous initial conditions in the energy landscape 
with a total energy barrier height of 5.569015. The diffusing particle is midway between the 
origin and destination sites while the bracing particles are displaced off the lattice line to 
accommodate the jump. Using molecular dynamics, we expect to relate this saddle point 
energy to a vacancy diffusion coefficient. 

The quintic potential can further be generalized by varying the distance at which the 
function is truncated. We set this cutoff distance to be v3. However, allowing this cutoff 
distance to vary introduces a larger class of potentials. A systematic study of the cutoff 
radius on the robustness of the ground states is necessary for experimental realization. The 
hard-core plus square shoulder potential has ground states that vary significantly depending 
on the relative lengthscale of the hard-core distance and the square-shoulder distance^^ 
It is expected that the ground states of the generalized quintic potential would be sensitive 
to the location of the minimum and the cutoff distance, though it is currently unknown 
how sensitive the ground-state phase diagram is to these parameters. Understanding this 
sensitivity is important to experimentalists who want a simple, robust potential. Developing 
an optimization procedure to make self-assembly more robust would be particularly useful. 
In addition, an extension to three-dimensions may provide additional fundamental insight. 

As mentioned earlier, inverse optimization techniques have been effective in develop- 
ing potentials for targeted material properties. We intend to develop a general and broad 
inverse optimization technique to target specific vacancy arrangements by accounting for 
and/or manipulating long-ranged vacancy- vacancy interactions. For example, one might de- 
velop an objective function whose variables include the strength, sign (attractive/repulsive), 
and angular dependence of vacancy-vacancy interactions. Alternatively, one might consider 
a two-component system system of a heavy particle and a light particle and apply an in- 
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verse optimization technique to this system. Using a broad family of potentials, one could 
then optimize over the available parameters to achieve a dilute concentration of effectively 
repulsive vacancies. 

VII. ACKNOWLEDGMENTS 

We thank Lawrence Cheuk for initial work on a related model. S.T. thanks the Institute 
for Advanced Study for its hospitality during his stay there. This work was supported by the 
Office of Basic Energy Sciences, U.S. Department of Energy, Grant DE-FG02-04-ER46108.. 



Corresponding author; Electronic address: torquato@electron.princeton.edu 



1 J. Sharma, R. Chhabra, A. Cheng, J. Brownell, Y. Liu and H. Yan, Science, 2009, 323, 112-116. 

2 A. P. Hynninen, J. H. J. Thijssen, E. C. M. Vermolen, M. Dijkstra and A. van Blaaderen, 
Nature Materials, 2007, 6, 202-205. 

3 Y. K. Cho, R. Wartena, S. M. Tobias and Y. M. Chiang, Adv. Func. Mat., 2007, 17, 379-389. 

4 S. Torquato, Soft Matter, 2009, 5, 1157-1173. 

5 H. Cohn and A. Kumar, Proc. Nat. Acad. Sci., 2009, 106, 9570-9575. 

6 W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal dispersions, Cambridge University 
Press: Cambridge, 1989. 

7 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2005, 95, 228301. 

8 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 73, 011406. 

9 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 74, 021404. 

10 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev, E, 2007, 75, 031403. 

11 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2008, 101, 85501. 

12 M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Phys. Chem. A, 2007, 111, 12816-12821. 

13 R. D. Batten, F. H. Stillinger and S. Torquato, J. App. Phys., 2008, 104, 033504. 

14 J. D. Joannopoulos, S. G. Johnson, J. N. Winn and R. D. Meade, Photonic crystals: molding 
the flow of light, Princeton University Press: Princeton, N.J., 2008. 

15 R. C. Agrawal and R. K. Gupta, J. Mat. Sci., 1999, 34, 1131-1162. 



25 



16 S. Torquato, Random Heterogeneous Materials: Microstucture and Macroscopic Properties, 
SpringerVerlag: New York, 2002. 

17 A. F. Andreev and L. M. Lifshitz, Sov. Phys. JETP, 1969, 29, 1107. 

18 E. Kim and M. H. W. Chan, Nature, 2004, 427, 225-227. 

19 N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders College: Philadelphia, Pa, 
1976. 

20 D. S. Fisher, B. I. Halperin and R. Morf, Phys. Rev. B, 1979, 20, 4692-4712. 

21 E. Cockayne and V. Elser, Phys. Rev. B, 1991, 43, 623-629. 

22 L. Candido, P. Phillips and D. M. Ceperley, Phys. Rev. Lett, 2001, 86, 492-495. 

23 L. Modesto, D. L. S. Junior, J. N. Teixeira Rabelo and L. Candido, Solid State Communications, 
2008, 145, 355-358. 

24 W. Lechner and C. Dellago, Soft Matter, 2009, 5, 646-659. 

25 W. Lechner and C. Dellago, Soft Matter, 2009, 5, 2752-2758. 

26 W. Lechner, E. Scholl-Paschinger and C. Dellago, J. Phys. Condens. Matter, 2008, 20, 404202. 

27 A. P. Hynninen, A. Z. Panagiotopoulos, M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. 
Chem. Phys., 2006, 125, 024505. 

28 M. Engel and H. R. Trebin, Phys. Rev. Lett, 2007, 98, 225505. 

29 M. Engel and H. R. Trebin, Zeitschrift fur Kristallographie, 2008, 223, 721-725. 

30 G. Doppelbauer, E. Bianchi and G. Kahl, J. Phys.: Condensed Matter, 2010, 22, 104105. 

31 J. Fornleitner and G. Kahl, Eur. Phys. Lett, 2008, 82, 18001. 

32 J. Fornleitner and G. Kahl, J. Phys: Condes. Matter, 2010, 22, 104118. 

33 M. A. Glaser, G. M. Grason, R. D. Kamien, A. Kosmrlj, C. D. Santangelo and P. Ziherl, Eur. 
Phys. Lett, 2007, 78, 46004:1-5. 

34 A. Quandt and M. P. Teter, Phys. Rev. B, 1999, 59, 8586-8592. 

35 M. Dzugutov, Phys. Rev. A, 1992, 46, 2984-2987. 

36 J. P. K. Doye, D. J. Wales and S. I. Simdyankin, Faraday Discussions, 2001, 118, 159-170. 

37 M. D. Haw, Phys. Rev. E, 2010, 81, 031402. 

38 D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, Inc. Orlando, 
FL, USA, 2001. 

39 B. J. Alder and T. E. Wainwright, Phys. Rev., 1962, 127, 359-361. 



26 



